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Abstract 

An asymptotic model based on a reductive perturbative expansion of 
the drift kinetic and the Maxwell equations is used to demonstrate that, 
near the instability threshold, the nonlinear dynamics of mirror modes in a 
magnetized plasma with anisotropic ion temperatures involves a subcritical 
bifurcation, leading to the formation of small-scale structures with ampli- 
tudes comparable with the ambient magnetic field. 

PACS: 52.35.Py, 52.25.Xz, 94.30.cj, 94.05.-a 



1. In regions of planetary magnetosheaths close to the magnetopause and 
in the solar wind as well, magnetic structures with a cigar form elongated along 
the direction of the ambient magnetic field are commonly observed (see e.g. [1, 
2J). According to recent observations [3], more than 60% of such structures are 
magnetic depressions (holes) associated with maxima of the density and pressure 
fluctuations. A typical depth of magnetic holes is about 20% of the mean magnetic 
field value and can sometimes achieve 50 %. The characteristic width of such 
structures is of the order of a few ion Larmor radii, and they display an aspect 
ratio of about 7-10. The origin of these structures is not fully understood but 
they are often viewed as associated with the nonlinear development of the mirror 
instability, a kinetic instability first predicted by Vedenov and Sagdeev pi] in 1957. 

The linear mirror instability has been extensively studied both analytically (see, 
e.g. [5j [6]), and by means of particle- in-cell (PIC) simulations [7J. This instability 
develops in a collisionless plasma, when the anisotropy of the ion temperature 
exceeds the threshold, 

TJT\\ - 1 = PZ 1 . (1) 

Here @± = 8np±/B 2 (similarly, = 8ttp^/B 2 ), where p± and pn are perpendicular 
and parallel plasma pressures respectively. Such conditions can be met under the 
effect of the plasma compression in front of the magnetopause [5]. As shown in 
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[SldniE], the instability is arrested at large k due to finite ion Larmor radius (FLR) 
effects. 

Mirror structures are also observed when the plasma is linearly stable [T2]. 
which may be viewed as the signature of a bistability regime. This property was 
also established in the framework of anisotropic magnetohydrodynamics, using an 
energetic argument [13J. The aim of the present paper is to demonstrate that 
the bistability of mirror structures results from a subcritical bifurcation. As well 
known, for such a bifurcation, non trivial stationary states below threshold are lin- 
early unstable, while above threshold, initially small-amplitude solutions undergo 
a sharp transition to a large-amplitude state, associated with a blowup behavior 
within an asymptotic theory. After reviewing the nonlinear theory of the mirror 
instability, briefly announced in [Hj, we demonstrate the subcritical character of 
the bifurcation in three steps: absence of small-amplitude stationary solution above 
threshold, existence of an unstable branch of solutions below threshold and blowup 
behavior for the initial value problem above threshold. 

The approach is based on a mixed hydrodynamic-kinetic description, assum- 
ing a weak nonlinear regime near threshold. Close to threshold, the unstable 
modes have wavevectors almost perpendicular to the ambient magnetic field B 
(k z /k± <C 1) with k±pi <C 1, so that the perturbations can be described using a 
long-wave approximation. The latter allows one to apply the drift kinetic equation 
(see, e.g., [HI [16]) to estimate the main nonlinear effects that correspond to a local 
shift of the instability threshold (1T1). All other nonlinearities connected, for exam- 
ple, with ion inertia are smaller. As the result, we obtain an asymptotic equation 
for the parallel magnetic field fluctuation, B z that displays a quadratic nonlinear- 
ity We show that this equation belongs to the generalized gradient type with a 
free energy that decreases in time, associated with the development of magnetic 
holes. This process has a self-similar blow-up behavior. This means that possible 
stabilization of the instability can only take place for amplitudes of order one, a 
regime that is beyond the framework of the present asymptotics. The present ap- 
proach contrasts with the quasi-linear theory [T7j that also assumes vicinity of the 
instability threshold but, being based on a random phase approximation, cannot 
predict the appearance of coherent structures. Phenomenological models based on 
the cooling of trapped particles were proposed to interpret the existence of deep 
magnetic holes [TSJ [T5] . These models do not however address the initial value 
problem in the mirror unstable regime. 

2. Consider for the sake of simplicity, a plasma with cold electrons. To describe 
the mirror instability in the long-wave limit it is enough to use the drift kinetic 
equation for ions ignoring parallel electric field En and transverse electric drift: 



In this approximation ions move along the magnetic field (b = B/5) due to the 
magnetic force /ib • VB where /x = vj_/ (2B) is the adiabatic invariant which plays 



df 
dt 



+ v\\b ■ Vf - fib ■ VB 




0. 



(2) 
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the role of a parameter in equation Both pressures pii and p± are given by 

Pll = mB J vjfd[idv\\d(p = m J v 2 fd 3 v, (3) 

p± = mB 2 J fj,fdfj,dv\\dip = -to y v 2 L fd 3 v. (4) 

Equation (J2J) with relations (j3J) and (jlj) are supplemented with the equation ex- 
pressing the balance of forces in a plane transverse to the local magnetic field 

4tt / n n B ■ VB i 

1+ m(^-fii)]-^r-} = °- ( 5 ) 



+ 



Here, consistently with the long-wave approximation, we neglect both the plasma 
inertia and the non-gyrotropic contributions to the pressure tensor. Furthermore, 
rijfc = 5ik — bibk denotes the projection operator in the plane transverse to the local 
magnetic field. In this equation, the first term describes the action of the magnetic 
and perpendicular pressures, the second term being responsible for magnetic lines 
elasticity. 

The equation governing the mirror dynamics is then obtained perturbatively 
by expanding Eqs. (jSJ) and (jSJ). In this approach, the ion pressure tensor elements 
are computed from the system (j2J), ([5]), near a bi-Maxwellian equilibrium state 
characterized by temperatures T± and Tji and a constant ambient magnetic field 
B taken along the z-direction. 

From Eq. (jSj) linearized about the background field B by writing B = B + B 
(B > B) with B ~ e - iLJ ' +ik ' r , we have 



(l) , B B Z _ k\ ( P±-P\\ \ B B Z 



Here k z and k± are the projections of the wave vector k, and pj^ is calculated from 
the linearized drift kinetic equation (121) : 

df^ + v df^_ dB z df® _ Q 
dt " dz dz dv\\ 

In Fourier space, this equation has the solution 

/W = -^-k,fH.. (7) 

UO — K 7 V\\ OV\\ 



The mirror instability is such that uj/k z Vth\\ = y2Tj|/TO. This means that the 
ions contributing to the resonance uj — k z v\\ = 0, correspond to the maximum of 
the ion distribution function. 
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After substituting ([7]) into the first order term for perpendicular pressure (pTJ 
and performing integration, we get 




iy/nuj (31 B Q B Z ^ 



\k z \v t h\\ P\\ 4vr 



The first term in (jHJ) is due to the difference between perpendicular and parallel 
pressures, while the second one accounts for the Landau pole. 

Equation (jSJ) together with ([6]) yield the growth rate for the mirror instability 
in the drift approximation where FLR corrections are neglected [I] 



7 = \kz\Vth\\ 



#l , 1 k 



2 



P\\ P± ki(3_ 



-X 



(9) 



where \ — 1 + (Ai — The instability takes place when /3±/f3\\ — 1 > /3J 1 

and, near threshold, develops in quasi-perpendicular directions, making the parallel 
magnetic perturbation dominant. 

As shown in Refs. [HI DUE], when the FLR corrections are relevant, the growth 
rate is modified into 



7 = \k z \vth\\ 



ki 3 



e - - —k 2 ±P 2 
k\ 4 X If 



(10) 



where e = /3±x 1 (P±/l3\\ — 1 — P^ 1 ) and the ion Larmor radius pi = v t ^/uj c 



is defined with the transverse thermal velocity Vth± = y2T±/m and the ion gy- 
rofrequency u C i = cBq / (mc) . This growth rate can be recovered by expanding the 
general expression given in [5] , in the limit of small transverse wavenumbers. It can 
also be obtained directly from the Vlasov- Maxwell (VM) equations in a long-wave 
limit which retains non gyrotropic contributions [20]. It is important to note that 
the expression (flOl) for 7 is consistent with the applicability condition u/k z <^ v t h\\, 
i.e. when the supercritical parameter \s\ 1. In this case the instability satura- 
tion happens at small k± oc ^fe due to FLR and for almost perpendicular direction 
in a small cone of angles, k z j k± oc ■sfe. As a result, the growth rate 7 oc e 2 , so 
that, when defining new stretched variables by 

k z = eK zP r\2/V3) X 1/2 , 

k ± = (2/V3)^K ± p- 1 x 1 / 2 , (11) 
7 = T(2/V3)e 2 Q (v^j.)- 1 (x/VA.) 3/2 , 



it takes the form 



T=\K Z \(1-K 2 /Kl-Kl). (12) 



Hence it is seen that, in the (K± — 0) plane (0 = K z /K±), the instability takes 
place inside the unit circle: 2 + K\ < 1. The maximum of T is obtained for 
K± = 1/2, = ± 1/2 and is equal to r max = 1/8. Outside the circle the growth 
rate becomes negative (in agreement with [TO]). 
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3. As it follows from in the linear regime, near the instability threshold, 
the fluctuations of perpendicular and magnetic pressures almost compensate each 
other (compare with (Q). Therefore, in the nonlinear stage of this instability, 
we can expect that the main nonlinear contributions come from the second order 
corrections to the total (perpendicular plus magnetic) pressure, i.e. 

(!) B B Z (2 ) B 2 d 2 B B Z 

P - + ^T + ^ + 8^ = - X A- ± —- (13) 

This result can be obtained rigorously by means of a multi-scale expansion based 
on the linear theory scalings ffTTj) . For this purpose, we introduce a slow time T 
and slow coordinates R in a way consistent with ffTTj) . and expand the magnetic 
field fluctuations as a powers series in e 1 ^ 2 : 

B z = eBP + 0(e 2 ), B ± = e 3 / 2 B ( l /2] + 0(e 5 / 2 ), (14) 

where B^/ 2 ) are assumed to be functions of R and T. Using these expressions 
it is easy to establish that quadratic nonlinear terms coming from the expansion 
of II in (J5j) as well as from the second term in the r.h.s. of Eq. (JT3l) are small 
in comparison with the quadratic term originating from the magnetic pressure in 
Eq. ([TBI . Thus, to get a nonlinear model for mirror dynamics it is enough to find 
Pj 2 ^. The expansion (|T4l) induces a corresponding expansion for the distribution 
function and for both pressures. Defining 

p^f 1 = Tim J v 2 L f ( - n ^v±dv±dvn, 

from (j4j) we have 

A 2) = {B?)/B f P f + 2(B^/B )^ + 

up to an additional contribution proportional to B^ that cancels out in the final 
equation due to the threshold condition. 

On the considered time scale, the effect of nonlinear Landau resonance is neg- 
ligible in the contribution to that can thus be estimated from the equation 

^ + (V/4,|)B»>-^^ = 0. 

For an equilibrium bi-Maxwellian distribution, we have 

fV = (2» 2 /vt hll )(BP) 2 f^ 

and thus 



P± 



R2 

( 2 ) _ (q AQl IQ i noi I n2^ °- 



(/? ± -4/? 2 //?|| + 3/3i//3j 



57T 



As a consequence, because of the vicinity to threshold we obtain 

P ? + fH 1 + ^)# >«. (15) 
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Then rewriting equation (1131) using the slow variables (Hip and rescaling the am- 
plitude 



we arrive at the equation 



dU ~ 
df = Kz 



(16) 



Here a = ±1, depending of the positive or negative sign of e, K z = —7~idz is a 
positive definite operator (whose Fourier transform is H is Hilbert transform 



1. 



f(Z') 



-dZ'. 



Hf(Z) = —VP 

71 J-oo Z 

As seen from the equation, its linear part reproduces the growth rate (fl2l . in 
particular, the third term in the r.h.s. accounts for the FLR effect. 

Equation ( fl6l) simplifies when the spatial variations are limited to a direction 
making a fixed angle with the ambient magnetic field. After a simple rescaling, 
one gets 



dU 
df 



U-3U 2 



;i7) 



where H is the coordinate along the direction of variation. This equation can be 
referred to as a "dissipative Korteveg-de Vries (KdV) equation" , since its station- 
ary solutions coincide with those of the usual KdV equation. The presence of 
the Hilbert transform in Eq. ffTTl) nevertheless leads to a dynamics significantly 
different from that described by soliton equations. 

4. Equation (fl6|) (and its ID reduction (fl7|) as well) possesses the remarkable 
property of being of the form 



dU 
df 



-K. 



6£ 
~5U' 



where 



- a -U>+ l -UA^dlU + l -{V 1 _Uf + U* 
-aN/2 + h/2 + I 2 /2 + h 



dR 



has the meaning of a free energy or a Lyapunov functional. This quantity can only 
decrease in time, since 



dF 
~dt 



6£dU_ 
6U~dt 



dR 



SF^ 5F 



K z —dR < 0. 



6U 5U 



(19) 



This derivative can only vanish at the stationary localized solutions, defined by the 
equation 

SF ( .^d 2 N 



6u = r _A± ^ 



+ A , \U-3U 2 



0. 



(20) 
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We now show that non-zero solutions of this equation do not exist above thresh- 
old (a = +1). For this aim, following Ref. [21], we establish relations between the 
integrals N, I\, I 2 and I3, using the fact that solutions of Eq. (1201) are stationary 
points of the functional F (i.e. 5F = 0). Multiplying Eq. (1201) by U and integrating 
over R gives the first relation 

aN -h-h- 34 = 0. 

Two other relations can be found if one makes the scaling transformations, Z — > aZ, 
H± — > &Rj_, under which the free energy (JTSJ) becomes a function of two scaling 
parameters a and b 

F ( a , b) = ~^ab 2 + ^a- 1 + |a + I 3 ab 2 . 

Due to the condition SF = 0, the first derivatives of F at a = b = 1 have to vanish: 

dF _ aN h.h 

d^-~~2~~ 2~ + 2" + /s " °' 
dF 

_ aiV + 2J 1 + 2/3 = 0. 

ob 

Hence, after simple algebra, one gets the three relations 

h + ^N = 0, J 3 = -2/ 1 , 72 = 3/!. 

For a — +1, the first relation can be satisfied only by the trivial solution U = 
0, because both integrals 7i and iV are positive definite. In other words, above 
threshold, nontrivial stationary solutions obeying the prescribed scalings do not 
exist. 

In contrast, below threshold, stationary localized solutions can exist. For these 
solutions, the free energy is positive and reduces to F s = N/2. Furthermore, I3 = 
/ U 3 d 3 R < 0. which means that the structures have the form of magnetic holes. As 
stationary points of the functional F, these solutions represent saddle points, since 
the corresponding determinant of second derivatives of F with respect to scaling 
parameters taken at these solutions, is negative (d aa FdbbF — {d^F) 2 = —2N 2 < 0). 
As a consequence, there exist directions in the eignefunction space, for which the 
free-energy perturbation is strictly negative, corresponding to linear instability of 
the associated stationary structure. This is one of the properties for subcritical 
bifurcations. 

For the ID model (fT7|) . the proof of instability of stationary solution Uq = 
— |sech 2 (S/2) (which coincides with the KDV soliton) is more complicated than 
in 3D. The corresponding free energy turns out to have a minimum relatively to 
the scaling parameter. Therefore, one needs to consider the linearized problem for 
perturbations W (U = Uq + W), which can be formulated as 

dW ^ 5F 



7 



where F = ^(W\L\W) is the quadratic part of the free energy and 

L = 1- — + 6U 

is the ID Schrodinger operator. 

It is easily seen that the operator L has one neutral (shift) mode ipx = dsUo 
(Ld^Uo = 0) associated with invariance by space translation, which has one node. 
Thus, according the oscillation theorem, L has one negative energy level with 
E = —5/4 < 0, corresponding to the ground state ip = sech 3 (S/2) (without 
nodes), which proves the instability of the stationary solution Uq with the growth 
rate equal l(ip \K\ip ) / (ip \ip ) > 0. 

As a consequence, starting from general initial conditions, the derivative dF/dt 
f|T9|) is almost always negative, except for unstable stationary points (zero measure) 
below threshold. In the nonlinear regime, negativeness of this derivative implies 
/ U 3 d 3 R < 0, which corresponds to the formation of magnetic holes. Moreover, 
this nonlinear term (in F) is responsible for collapse, i.e. formation of singularity 
in a finite time. 

5. In order to characterize the nature of the singularity of Eq. (TP71) . it is 
convenient to introduce the similarity variables £ = (T — T) _1 / 3 S, r = — log(T ( 



o 



T), and to look for a solution in the form U = (T — T) 2 ^ 3 g(£, r), where g(£, r) 
satisfies the equation 



dg 2 £dg v 



d 2 g 2 

w~ 39 



As time T approaches T (r — > oo), the last term in this equation becomes neg- 
ligibly small and simultaneously d T g —>■ so that asymptotically the equation 
transforms into 

( \d 2 9 2 

TP? ^9 



s 9+ 3d-r K ' 



(21) 



For the free energy this means that close to To the first term ~ iV turns out to be 
much smaller in comparison with all other contributions, in particular with J £/ 3 cG. 

At large |£|, that corresponds to the limit T — > T , the asymptotic solution g 
of Eq. (HU) obeys 

d£ 

where C = J JZc g 2 (t')d£' > 0, and has the form g= C^ 2 log |£/£ |- For U, it 
gives the asymptotic solution 

C 

U asyrnp = ^ log |E/Eo(t)| 

with S (t) = (Tq — T) 1 / 3 ^, that, as T — > T , has an almost time independent 

1 /3 

tail. For |S| < (T — T) |£ | the solution is negative and becomes singular as H 
approaches the origin. 



S 



6. Asymptotically self-similar solutions can also be constructed in three dimen- 
sions, when rescaling the longitudinal coordinate by (T — T) 1 / 2 , the transverse ones 
by (To — T) 1//4 and the amplitude of the solution by (To — T) -1 / 2 . Existence of a 
finite time singularity for the initial value problem can be established for initial 
conditions for which the functional F is negative, when the term involving a can 
be neglected, an approximation consistent with the dynamics: 

F - T lim = | + | + 7 3 . (22) 



To prove this statement, consider the operator K~ l , (inverse of the operator K 
which is defined on functions obeying J U(Z,H±)dZ = 0, a condition consistent 
with Eq. (fl6l) . Then the time derivative of Fn m can be rewritten through the 
operator K~ x as follows, 



dF 



lim 



dT 



J UrK^UrdR < 0. (23) 



Consider now the positive definite quantity N = J UK Z 1 UdH > 0, whose dynamics 
is determined by the equation 

dN 

— = -2 {h + I 2 + 3/ 3 ) = -6F lim + h + I 2 . (24) 

Let Ti; m be negative^ initially, then at T > the r.h.s. of ( |24l) will be positive, and, 
as a consequence, N will be a growing function of time. 

Introduce now the new quantity S = —Fu m /N which is positive definite if 
Fum\T=o < 0. The time derivative of S is then defined by means of Eqs. (1231) and 
(El: _ 

lIS FhmNT + 1 { U T K?U T dBL (25) 



dT N 2 N 

The second term in the r.h.s. of this equation can be estimated using the Cauchy- 
Bunyakowsky inequality: 



dN r / r \ 1 / 2 

— = 2 J UK; l U T d~R < 2N 1/2 (J U T K- x U T d^\ 



that gives 

J U T K; l U T dK > iV|/(4iV). 

Substituting the obtained estimate into Eq. (1251) and taking into account definition 
(T22|) for F\ im and Eq. (1241 as well, we arrive at the differential inequality for S 
(compare with [22]): 



dS N T 
dT ~W 



Nt_ f 
4 lim 



> 15 S 2 . 



Integrating this first-order differential inequality yields 



S > 1 ■ (26) 
15 (-To — T) 
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Here the collapse time T = (15 So) -1 is expressed in terms of the initial value 
S\t=o = So- It is interesting to mention that the time behavior of S given by the 
estimate (1261) coincides with that given by the self-similar asymptotics. 

7. In this letter, we have presented an asymptotic description of the nonlinear 
dynamics of mirror modes near the instability threshold. Below threshold, we have 
demonstrated the existence of unstable stationary solutions. Differently, above 
threshold, no stationary solution consistent with the prescribed small-amplitude, 
long-wavelength scaling can exist. For small-amplitude initial conditions, the time 
evolution predicted by the asymptotic equation ([16]) leads to a finite-time singu- 
larity. These properties are based on the fact that this equation belongs to a class 
of generalized gradient systems for which it is possible to introduce a free energy 
or a Lyapunov functional that decreases in time. In ID, this model can be referred 
to as the dissipative KDV equation. The difference with the usual KDV equation 
is connected with the change of the symplectic operator d/dZ for the KDV case 
to the positive definite operator K z = —Hd/dZ for (ITTj) . This leads to significant 
changes in the system dynamics. 

The singularity formation as well as the existence of unstable stationary struc- 
tures below the mirror instability threshold obtained with the asymptotic model, 
can be viewed as features of a subcritical bifurcation towards a large-amplitude 
state that cannot be described in the framework of the present analysis. Such 
an evolution should indeed involve saturation mechanisms that become relevant 
when the perturbation amplitudes become comparable with the ambient field. We 
note that recent numerical simulations [231 ED] of the VM equations, using either 
particle-in-cell codes or an Eulerian description, display the formation of mag- 
netic humps above threshold together with a phenomenon of bistability, associated 
with the existence of stable large-amplitude magnetic holes both below and above 
threshold. Due to numerical constraints, these simulations are however performed 
in regimes that are not sufficiently close to threshold, for the present theory to be 
applicable. 

An important question concerns the relation between the present theory of 
structure formation and the quasi-linear effects that could compete near threshold. 
An early-time quasi-linear regime could for example modify the onset of coherent 
structures and, on the other hand, the development of such structures, can also 
affect the quasi-linear dynamics. 

8. This work was performed in the framework of ISSI team "The effect of 
ULF turbulence and flow chaotization on plasma energy and mass transfers at the 
magnetopause" . The work of EK was supported by RFBR (grant no. 06-01-00665) 
and by the French Ministere de l'Enseignement Superieur et de la Recherche during 
his visit at the Observatoire de la Cote d'Azur, that of TP and PLS by "Programme 
National Soleil Terre" of CNRS. 
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